
data{
  int<lower=1>K; //number of conefficients
  int<lower=1>N; //number of observations
  matrix[N,K] x;
  int<lower=1> C; #n counties
  int<lower=1> P; #n parties
  int<lower=1, upper=C> county[N];
  int<lower=1, upper=P> party[N];
  vector<lower=0, upper=100>[N]y;
}

transformed data{
matrix[N,K] Qx=qr_Q(x)[,1:K]*sqrt(N-1);
matrix[K,K] Rx=qr_R(x)[1:K,]/sqrt(N-1);
matrix[K,K] Rx_inv=inverse(Rx);
}



parameters{
  vector[K] beta_tilde;
  vector[C] ranef_c;
  vector[P] ranef_p;
  real intercept;
  real<lower=0, upper=100>sigma_y;
  real<lower=0, upper=100>sigma_b;
  real<lower=0, upper=100>sigma_c;
  real<lower=0,upper=100>sigma_p;
}

transformed parameters{
  vector[N] eta;
  vector[N] y_hat;
  vector[K] beta=Rx_inv*beta_tilde;
  eta=Qx*beta_tilde;


  for(n in 1:N){
    y_hat[n]=intercept + eta[n];
  }

}

model{
    sigma_c~inv_gamma(5,5);
    sigma_p~inv_gamma(5,5);
    intercept~normal(0,1);
    to_vector(ranef_c)~normal(0,sigma_c);
    to_vector(ranef_p)~normal(0,sigma_p);

  for(n in 1:N){
    y[n]~normal(y_hat[n], sigma_y);
  }


}

generated quantities{
  vector[N] y_rep;
  vector[N] log_lik;
  vector[K] fd_beta;

  for(n in 1:N){
    log_lik[n]=normal_lpdf(y[n]|y_hat[n], sigma_y)  ;
  }

  for(n in 1:N){
    y_rep[n]=normal_rng(y_hat[n], sigma_y);
  }

  for(k in 1:K){
    if(k<6)
      fd_beta[k]=(beta[k]*1)-(beta[k]*0);
    else
      fd_beta[k]=(beta[k]*(mean(x[,k]+sd(x[,k]))))-(beta[k]*(mean(x[,k]-sd(x[,k]))));
      }
}
